home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
FishMarket 1.0
/
FishMarket v1.0.iso
/
fishies
/
201-225
/
disk_214
/
mandelvroom
/
src
/
mandffp.c
< prev
next >
Wrap
C/C++ Source or Header
|
1992-05-06
|
9KB
|
426 lines
/*
* MandelVroom 2.0
*
* (c) Copyright 1987,1989 Kevin L. Clague, San Jose, CA
*
* All rights reserved.
*
* Permission is hereby granted to distribute this program's source
* executable, and documentation for non-comercial purposes, so long as the
* copyright notices are not removed from the sources, executable or
* documentation. This program may not be distributed for a profit without
* the express written consent of the author Kevin L. Clague.
*
* This program is not in the public domain.
*
* Fred Fish is expressly granted permission to distribute this program's
* source and executable as part of the "Fred Fish freely redistributable
* Amiga software library."
*
* Permission is expressly granted for this program and it's source to be
* distributed as part of the Amicus Amiga software disks, and the
* First Amiga User Group's Hot Mix disks.
*
* contents: this file contains the routines to calculate Mandelbrot and
* Julia pictures in Amiga Fast Floating Point.
*/
#include "mandp.h"
extern SHORT MaxOrbit;
static LONG MaxIter;
/*
* Fast Floating Point Mandelbrot Generator
*/
MandelbrotFFP( Pict, StartX, StartY, GapX, GapY )
/*
* These parameters are IEEE 32 bit floating point numbers
*/
struct Picture *Pict;
LONG StartX, StartY, GapX, GapY;
{
register int i, j, k;
register SHORT *CountPtr;
register float curx, cury;
float startx, starty, gapx, gapy;
float SPFieee();
UBYTE MandFlag;
/* Here we convert the IEEE parameters into FFP format */
startx = SPFieee( StartX );
starty = SPFieee( StartY );
gapx = SPFieee( GapX );
gapy = SPFieee( GapY );
MaxIter = Pict->MaxIteration;
if (Pict->Flags & NO_RAM_GENERATE)
CountPtr = Pict->Counts;
else
CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;
/* start in the upper left hand corner */
cury = starty;
/*
* for each row
*/
for (i = Pict->CurLine; i < Pict->CountY; i++) {
curx = startx;
MandFlag = 0;
if ( Pict->Flags & NO_RAM_GENERATE )
CountPtr = Pict->Counts;
/*
* for each collumn
*/
for (j = 0; j < Pict->CountX; j++) {
/* if we've hit a Mandelbrot point, then use the convergence detector,
* to reduce compute time
*/
if (*CountPtr == 0) {
if ( MandFlag ) {
k = FFP_Trace_Height( 0.0, 0.0, curx, cury );
} else {
k = FFP_Height( 0.0, 0.0, curx, cury ); /* Normal calculator */
}
MandFlag = k == Pict->MaxIteration;
*CountPtr = k; /* save it in ram for recoloring */
}
CountPtr++;
curx += gapx;
ChildPause( Pict );
}
cury += gapy;
CheckEOL( Pict );
}
} /* Mandelbrot */
/***************************************************************************
*
* Julia Curve Generator that uses Amiga Fast Floating Point
*
**************************************************************************/
/*
* Fast Floating Point Julia Generator
*/
JuliaFFP( Pict, JuliaX, JuliaY, StartX, StartY, GapX, GapY )
/*
* These parameters are IEEE 32 bit floating point numbers
*/
register struct Picture *Pict;
LONG JuliaX, JuliaY;
LONG StartX, StartY;
LONG GapX, GapY;
{
register int i, j, k;
register float curx, cury;
register SHORT *CountPtr;
float juliax, juliay, gapx, gapy,startx;
float SPFieee();
UBYTE ConvFlag;
MaxIter = Pict->MaxIteration;
if (Pict->Flags & NO_RAM_GENERATE)
CountPtr = Pict->Counts;
else
CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;
/* Here we convert the IEEE parameters into FFP format */
juliax = SPFieee( JuliaX );
juliay = SPFieee( JuliaY );
gapx = SPFieee( GapX );
gapy = SPFieee( GapY );
/* start in the upper left hand corner */
startx = SPFieee(StartX);
cury = SPFieee(StartY);
/* move down to next not generated line */
cury += Pict->CurLine*gapy;
/*
* for each row
*/
for (i = Pict->CurLine; i < Pict->CountY; i++) {
curx = startx;
ConvFlag = 0;
if ( Pict->Flags & NO_RAM_GENERATE )
CountPtr = Pict->Counts;
/*
* for each collumn
*/
for (j = 0; j < Pict->CountX; j++) {
/* if we've hit a Mandelbrot point, then use the convergence detector,
* to reduce compute time
*/
if (*CountPtr == 0) {
if (ConvFlag) {
/* try the convergence method */
k = FFP_Trace_Height( curx, cury, juliax, juliay);
} else {
k = FFP_Height( curx, cury, juliax, juliay); /* Normal calc */
}
ConvFlag = k == Pict->MaxIteration;
*CountPtr = k; /* save it in ram for recoloring */
}
CountPtr++;
curx += gapx;
ChildPause( Pict );
}
cury += gapy;
CheckEOL( Pict );
}
} /* JuliaFFP */
/*
* This function calculated the 'potential' of a given location
* in the complex plane.
*/
int
FFP_Height( posx, posy, juliax, juliay )
float posx; /* Location in screen coordinate space */
float posy; /* Location in screen coordinate space */
float juliax; /* Real part of location on complex plane */
float juliay; /* Imaginary part of location on complex plane */
{
register float cura, curb;
register float cura2, curb2;
register LONG k;
#ifdef CHECK_TASK_STACK
CheckStack();
#endif
cura = cura2 = posx;
curb = curb2 = posy;
cura2 *= cura2;
curb2 *= curb2;
for (k = 0; k < MaxIter; k ++ ) {
curb *= cura; /* b = 2 * a * b + juliay */
curb += curb + juliay;
cura = cura2 - curb2 + juliax; /* a = a2 - b2 + juliax */
cura2 = cura * cura; /* a2 = a * a */
curb2 = curb * curb; /* b2 = b * b */
if (cura2+curb2 >= 16.0) /* quit if a2+b2 is too big */
return( k );
}
return( k );
}
/*
* this is also a Mandelbrot potential calculator, that is tailored to
* identify points that are 'in' the Mandelbrot set. Points outside the
* Mandelbrot set diverge. Points inside converge.
* This code uses a trace table mechanism to detect convergence.
*/
int
FFP_Trace_Height( posx, posy, curx, cury )
float posx; /* Real part of location on complex plane */
float posy; /* Imaginary part of location on complex plane */
float curx; /* Real part of location on complex plane */
float cury; /* Imaginary part of location on complex plane */
{
LONG k, l;
LONG TraceSize = 32;
float olda[32], oldb[32]; /* Convergence trace table */
register float cura, curb, cura2, curb2;
register float *RealTrace, *ImagTrace;
#ifdef CHECK_TASK_STACK
CheckStack();
#endif
cura = cura2 = posx;
curb = curb2 = posy;
cura2 *= cura2;
curb2 *= curb2;
for (k = 0; k < MaxIter; k += TraceSize) {
RealTrace = olda;
ImagTrace = oldb;
for (l = 0; l < TraceSize; l++) {
*(RealTrace++) = cura;
*(ImagTrace++) = curb;
curb *= cura;
curb += curb + cury;
cura = cura2 - curb2 + curx;
cura2 = cura * cura;
curb2 = curb * curb;
if (cura2+curb2 >= 16.0)
return( k + l );
}
/* Scope out trace table for convergence */
RealTrace = olda;
ImagTrace = oldb;
for (l = 0; l < TraceSize; l++)
if (cura == *(RealTrace++) && curb == *(ImagTrace++)) {
k += MaxIter;
return( MaxIter );
}
}
return( MaxIter );
}
DrawOrbitFFP( Pict, creal, cimag, zreal, zimag )
register struct Picture *Pict;
LONG zreal, zimag, creal, cimag;
{
register struct RastPort *Rp;
register float cura, curb;
float cura2, curb2;
float Creal, Cimag;
float x_scale, y_scale;
int x_center, y_center;
int width, height;
int x, y;
register int k;
struct Window *Window;
Window = OrbitWind;
Rp = Window->RPort;
width = (Window->Width-Pict->LeftMarg-Pict->RightMarg);
height = (Window->Height-Pict->TopMarg-Pict->BotMarg);
x_center = width/2 + Pict->LeftMarg;
y_center = height/2 + Pict->TopMarg;
y_scale = x_scale = (float) height / 2.0;
/*x_scale *= AspectRatio( Pict );*/
if ( Pict->pNode.ln_Type == MANDPICT ) {
Creal = SPFieee(creal);
Cimag = SPFieee(cimag);
cura = cura2 = SPFieee(zreal);
curb = curb2 = SPFieee(zimag);
} else {
Creal = SPFieee(zreal);
Cimag = SPFieee(zimag);
cura = cura2 = SPFieee(creal);
curb = curb2 = SPFieee(cimag);
}
cura2 *= cura2;
curb2 *= curb2;
SetAPen( Rp, 0 );
RectFill( Rp, Pict->LeftMarg, Pict->TopMarg, width+Pict->LeftMarg-1, height+Pict->TopMarg-1);
SetAPen( Rp, HIGHLIGHTPEN);
for (k = 0; k < MaxOrbit; k++ ) {
curb *= cura;
curb += curb + Cimag;
cura = cura2 - curb2 + Creal;
cura2 = cura * cura;
curb2 = curb * curb;
if (cura2+curb2 >= 16.0)
return( k );
/* map real and imaginary parts into window coordinates */
x = x_center + (int)(x_scale*cura);
y = y_center + (int)(y_scale*curb);
if ( x >= Pict->LeftMarg && x < Window->Width-Pict->RightMarg &&
y >= Pict->TopMarg && y < Window->Height-Pict->BotMarg ) {
/* plot pixel location */
WritePixel( Rp, x, y);
}
}
return( k );
}